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A  NUMERICAL  SOLUTION  TO  AN  ABLATION  PROBLEM  WITH 
POSSIBLE  LASER  APPLICATIONS 

ABSTRACT 

The  problem  of  an  ablating  solid  of  limited  extent  under  an  applied 
non-uniform  heat  input  is  solved  numerically,  and  the  solution  is  com¬ 
pared  in  certain  limiting  cases  with  results  from  the  exact  analytical 
equations.  Further,  the  usefulness  of  this  model  in  predicting  the 
depth  of  hole  made  in  materials  by  a  laser  beam  is  explored. 
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1.  INTRODUCTION 


The  problem  of  an  ablating*  solid  of  limited  extent  under  an 
applied  non-uniform  heat  input  is  solved  numerically,  and  the  solution 
is  compared  in  certain  limiting  cases  with  results  from  the  exact 
analytical  equations.  Further,  the  usefulness  of  this  model  in  predict¬ 
ing  the  depth  of  hole  made  in  materials  by  a  laser  beam  is  explored. 

Problems  similar  to  the  one  considered  have  been  formulated  by 

1_8** 

numerous  investigators,  and  some  particular  solutions  exist.  A 

numerical  solution  for  this  problem  has  been  given  by  Landau'*'  for  the 

special  case  of  a  constant  heat  input  to  a  semi -infinite  solid.  More 

k 

recently  Ready  reporting  the  effects  of  the  absorption  of  laser 
radiation  on  metals  introduced  results  from  the  ablation  model  for  the 
case  of  a  semi -infinite  solid  under  a  non-uniform  heat  input;  but,  he 
gave  no  details  of  the  solution.  In  view  of  the  very  restrictive 
(specialized)  treatment  which  this  problem  has  received  in  the  litera¬ 
ture  and  because  of  its  possible  application  to  certain  laser  effects, 
we  find  sufficient  justification  in  its  reconsideration. 

2.  STATEMENT  OF  THE  PROBLEM 

A  slab  of  thickness  B  initially  at  temperature  TQg(x)  is  heated 
on  one  end  with  energy  Hf(t)  over  an  area  with  dimensions  large  enough 
in  relation  to  the  final  depth  of  heat  penetration  that  heat  flow 
becomes  mathematically  one -dimensional.  The  other  end  of  the  slab  is 
insulated. 


Ablation  is  generally  associated  with  vaporization  and  sublimation , 
but  the  meaning  has  been  extended  to  include  liquefaction  in  which 
the  liquid  is  -removed  as  it  is  formed. 

** 

Superscript  numbers  denote  references  which  may  be  found  on  page  27. 
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Given  a  pulse  of  sufficient  power  density  and  a  long  enough  time, 

* 

the  surface  of  the  slab  will  reach  the  phase -change  temperature.  The 

material  which  has  been  transformed  is  removed  immediately  and  the  sur¬ 
face  recedes  into  the  solid.  We  are  interested  primarily  in  determining 
the  thickness  of  material  transformed  with  time. 

We  make  the  following  definitions: 

c  =  specific  heat  of  material  at  constant  pressure,  (cal/g°C) 
p  =  density  of  material,  (g/cm^) 

K  =  thermal  conductivity  of  material,  (cal/cm  °C  sec) 

O 

K  =  thermal  diffusivity  of  material,  equal  to  K/cp,  (cm^/sec) 

L  =  latent  heat  of  transformation  of  material,  (cal/g) 

T  g(x)  =  initial  temperature  distribution  of  material,  (°C);  for 
°  g(x)  =  1  the  temperature  is  initially  uniform  throughout 

the  material. 

T  =  transformation  temperature  of  material,  (°C) 

y(t)  =  thickness  of  material  transformed,  (cm) 

Hf(t)  =  heat  input  as  a  function  of  time  where  H  is  a  given  constant, 
o 

(cal/cm^sec) ;  f(t)  =  1  for  a  uniform  heat  input 


B  =  initial  thickness  of  material,  (cm) 

The  following  equations  describe  the  process : 

Tt  =  KTxx  y(t)  -  x  <  B,  t  >  0  ,  (2.1) 

T(x,0)  =  T  g(x)  <  Tt  0  <:  x  <:  B  ,  (2.2) 

O  .L» 

Tx(B,t)  =  0  ,  (2.3) 

Hf(t)  =  -  KTx  y(t)  =  0  ,  (2.4) 


Whether  the  phase-change  temperature  is  that  for  melting  or  vaporiza¬ 
tion  will  depend  on  the  incident  power  density.  At  low  densities  only 
liquid  will  appear ;  whereas ,  at  high  densities  the  liquid  formed  prior 
to  vaporization  can  he  ignored  because  the  latent  heat  of  melting  is 
much  smaller  than  that  of  vaporization. 
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(2.5) 


Hf(t)  =  -  KTx  +  pL  M|1 

M^l  s  0  for  T(y(t),t)  =  TL  ,  (2.6) 

at  -L* 

Mil  =  0  for  T(y(t),t)  <  T.  .  (2.7) 

at  -Lj 

Equation  (2.1)  may  be  recognized  as  the  general  parabolic  relation 
of  heat  flow.  The  initial  temperature  distribution  of  the  slab  is  given 
by  (2.2),  while  (2.3)  expresses  the  condition  of  no  heat  losses  from  the 
backface.  Equation  (2.5)  is  a  statement  relating  the  heat  input  to  the 
rate  of  heat  flow  into  the  solid  plus  the  rate  of  heat  absorbed  in  the 
transformation.  Tx  is  evaluated  at  the  position  of  the  boundary  y(t). 
Before  the  phase -change  temperature  is  reached  there  is  no  boundary 
motion,  (2.7),  so  that  (2.5)  reduces  to  (2.4).  Equation  (2.6)  states 
that  when  the  surface  reaches  the  temperature  T^  there  may  or  may  not 
occur  boundary  motion;  this  will  depend  on  whether  the  energy  required 
for  the  phase-change  is  or  is  not  absorbed  at  the  surface. 

The  equations  can  be  manipulated  more  readily  if  the  following 
transformations  are  made: 

,  T(x,t) 

v(x,t)  =  V  5 

L 

T  =  Kt  ) 

N  =  gp-  > 

CXL 

H 

p  =  ktl  • 

Equations  (2.l)  through  (2.7)  then  become 

v  =  v  y(T)  5?x<B>  t>0  ,  (2.8) 

T  xx 

v(x,0)  =  v^g(x)  <1  O^x^B*  (2.9) 
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V  (B,t)  =0  ,  (2.10) 

X 

pf (t )  =  -  v  y(r )  =  0  ,  (2.11) 

pf(r)  =  -  vx+  N&£ll  ,  (2.12) 

^  0  for  v(y(T),T)  =  1  ,  (2.15) 

=  0  for  v(y(T),T)  <  1  •  (2.14) 


5.  FINITE  DIFFERENCES 

We  denote  v.  .  to  be  the  value  of  v  at  the  point  i  Ax  at  time  j  At 
j  <3 

and  v.  .  ,  the  corresponding  value  at  time  ( j+l)  At.  We  further  define 
J"*'— 

according  to  practice  h  =  Ax,  k  =  At  and  r  =  — p  .  If  the  Crank - 

6  h 

Nicholson  differencing  scheme^  is  selected,  the  time  derivative  can  be 
approximated  by 


=  vi,j+i :  vi,j 


— 

0  1 


2  n 


(5.1) 


and  the  second  order  space  derivative  by 


xx 


~  2h2  K  vi+1, 


+  v  ^  -  2  f  v,  .  ,  +  v .  .) 

j+l  i+l,j  V  i,J+l  i,J./ 


+  ^  Vi-1, j+l  +  Vi- 


(5.2) 


where  0  £  J  is  the  truncation  error.  Substituting  (3.1)  and  (3«2)  in 

V 

(2.8),  introducing  r  =  ^  and  rearranging  we  obtain 

h 


2  Vi-1,  j+l  +  (1+r)vi^  j+i  "  2  Vi+1,  j+l  = 


!Vi-l,j+  (l_r)vi,j  +  !  Vi+1,, 


(3.3) 
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Equation  (3-5)  does  not  hold  at  the  mathematical  surface  i  =  0  since 

no  account  has  been  taken  of  the  energy  input  to  the  material.  This  is 

dealt  with  by  expanding  v.  .  in  a  Taylor  series  about  the  point  (0,  j+l) 

1  >  J 

in  both  the  space  and  the  time  directions.  Thus 


'l,j+l  V0,j+1  '  LdxJO,j+l 

(Ax)5  1  • 


r§£Ui  ^ + ^  ^ 


+  o 


V  =  V 

0,j  0 


r|vl  At  +  T  (At  ) 

,  j+1  LStJo^j+I  L 


2  1 


(3-4) 


(3-5) 


then  from  (3-4), 


is  replaced  by  [Sv/Bt]q  in  accordance  with 
(5.5),  and  (2.1l)  the  resulting  equation  becomes 


(2.8), 


(1+2r>T0,^i  -  -  vo,o  +  2rhpf+i 


(3.6) 


Equation  (3*6)  applies  up  to  the  time  that  the  phase-change  temperature 
is  reached,  i.e.,  so  long  asy(T)  =  0. 

Beyond  this  point  account  must  be  taken  of  the  fact  that  the 

boundary  no  longer  remains  stationary.  To  avoid  ambiguity  in  the  use  of 

the  space  coordinate  i,  a  new  coordinate  w  is  defined  such  that  the 

boundary  will  always  lie  at  w  or  between  w  and  w+  1.  At  the  original 

surface  w  =  0.  The  index  w,  an  integer,  depends  upon  time,  hence  we 

denote  this  dependence  by  w^.  Initially,  that  is,  at  j  =  0,  w^  —  0, 

however,  at  later  times  0  <  v.  <  w^n  and  for  each  j  the  boundary  lies 

between  w  and  w  , .  During  boundary  motion  the  temperature  of  the 
j  j+1 

surface  according  to  (2.13)  is  a  known  constant  so  that  the  first 

unknown  temperature  lies  at  the  grid  point  w  +  1  and  is  designated 

-,r  .  We  exnand  about  this  point  in  both  the  time  and  space 

"w+l,  j+1’ 

directions, 


r»vi 


vw+s,j+l  Vw+l,j+l  LdxJw+1,  j+1 

.2 


(l-s)  Ax 


(1/2)  m 

dx 


w+1,  j+1 


(l-s)2Ax^  +  0 


[(Ax)5  ]  , 


(3-7) 


ll 


(3.8) 


vw+l,j  vw+l,j+l 


vw+2,j+l  vwfl 


[IrL+i.j+i^  +  0  [(At)2  J  ’ 
.  ,  +  ri^i  ,  .  ,  ax 

,J+X  LOXJW+i.,  J-KL 


(3.9) 


and  where  0  £  s  £  1.  But  at  the  surface  we  have  v  .  ,  =  1.  Further, 

w+s,j+l  > 

if  (5*7)  and  (3.9)  are  added  after  (3*9)  has  been  multiplied  through  by 
(l-s)  and  (2.8)  and  (3*8)  are  substituted  in  that  order  in  the  resulting 
equation,  we  obtain 


(2-s>])  (2r  +  Vi, 

C1'8 0+l3  C2"8 j  +  2r 


j+1  ■  2r^.1‘,J+l'*  Vw+2,JH 


(5.10) 


For  clarity  the  subscript  j+1  has  been  attached  to  s  to  show  the 
exact  time  at  which  this  quantity  is  evaluated. 

Next  the  position  of  the  boundary  y(T)  must  be  determined.  From 
the  preceding  discussion  it  is  clear  that  y  may  be  expressed  in  terms 
of  an  integer  w  and  a  fraction  s  by 


yo = (wj + s? h  - 


and 


(3.11a) 


=  (vi +  vi)  h  • 


(3.11b) 


Subtracting  (3.11a)  from  (3.11h)  we  obtain 


yju  ’  ^ =  h  [Cvx  •  wi) +  0>i  •  sj)l •  (5.12) 
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fi-iru-'p  v.  w  and  s  are  zero  at  the  start  of  the  transformation,  there 
J  J  j 

is  no  difficulty  in  keeping  track  of  the  boundary  for  subsequent  times 
provided  Ay  can  be  calculated.  This  may  be  done  from  (2.12)  which  has 
the  finite  difference  form 

k^  =  [Vx  +  pf j+ll  /N  ’  (5.13) 

where  v  is  the  approximation  of  the  space  derivative  v  at  the  surface 

X 

and  may  be  evaluated  by  a  three  point  Gregory -Newt on  forward  difference 
approximation  as 

vx  =  L-  (v5'2sj+l.')  vw+l,  j+1 

+  4C2  -  8j+l)  V^2,jtl  (3.14) 

"  2sj+l^  vw+3,  jt-lJ  ^2h 

Equation  (3. 13)  can  be  arranged  more  conveniently  in  terms  of  r  giving 

T  -  f  [\  *  jd  •  (5-15) 


One  is  confronted  in  (3*15)  with  the  evaluation  of  s  from  temperatures 

j+i 

which  are  not  known  and  which  themselves  cannot  be  calculated  without 
knowing  s  .  ..  .  To  circumvent  this  problem  an  iterative  procedure  must  be 

j+1 

used. 

Finally,  we  note  that  at  the  back  surface  (3*3)  must  be  modified 

to  take  into  account  the  boundary  condition  (2.10).  Formally,  this  is 

* 

done  by  extending  the  mesh  one  unit  past  the  back  face,  where  the 
latter  is  indexed  by  i  =  M  =  B/Ax.  The  existence  of  a  zero  temperature 
gradient  requires  that  vw  n  =  vMi1  which  reduces  (3*3)  at  the  back  face 

1VXTJ_ 

to 


■"M-l,  J.1  +  (1  +  =  (1  •  r)V,M,j  +  "M-l.y  (5-l6> 


Creating  a  false  boundary. 
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4.  CALCULATING  PROCEDURE 


4 ■ 1  General 

The  difference  equations  of  Section  3  from  which  the  temperatures 
are  computed  have  the  tridiagonal  appearance  (for  each  mesh  row) 


B0V0  +  Vl  =  d0 


Alv0  +  Blvl  +  C1V2  =  dl 


0vQ  +  A2v1  +  B?v2  +  =  d2 


(4.1) 


+  BM-1VM-1  +  CM-1VM  “  Vl 


amvm-i  +  bmvm  Si  ’ 


with  zeros  everywhere  except  on  the  main  diagonal  and  on  the  two 
diagonals  parallel  to  it  on  either  side.'7  In  (4.l)  all  known 
quantities  have  been  lumped  in  d  thereby  eliminating  the  necessity  of 
the  further  use  of  j.  Through  a  Gaussian  elimination  process  this 
system  can  be  solved  explicity  for  the  unknown  v’s.  The  method  is 
credited  to  L.  H.  Thomas but  received  wide  attention  in  published 
form  in  an  article  by  Bruce,  Peaceman,  Rachford,  and  Rice.dd  Expressed 
in  a  less  cumbersome  waydd  the  system  (4.1)  becomes 


A.  v 

l 


i-1 


V 


C0V1  =  d0  ’ 

(a) 

iVi  *  di  1  s  s  M  - 1’ 

( 4 . 2 )  ( b ) 

bmvm  =  ^  ’ 

(c) 
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kst 


where  the  unknown  temperatures  corresponding  to  the  (j+l)  step  are 

given  by 


with 


VM  ~  % 

Vi  =  h.  '  Vi+l 


ona-i, 


4°“  Bo 


*i  '  VVi-1 


(**.3) 


*  _  10 
°"  Bo 


bi  "  B. 


'1  ! 


AJj 


ih-i 

I 

I 


1  S  i  S  M  ■  1  . 


For  convenience  the  calculations  are  carried  out  in  two  parts  - 

i 

Part  I  covering  up  to  the  time/  the  phase -change  temperature  is  reached 
and  Part  II  beyond  this  point.  In  Part  I  (4.2)  represents  the  system 
(3-6),  (3.3)  and  (3.16);  in  Part  II  (4.2)  represents  the  system  (3.10), 
(3.3)  and  (3.16)  with  (3.15)  ind  (3*12)  being  used  in  locating  the 
boundary.  Also  in  Part  II,  (4.2a)  is  more  appropriately  indexed  as 


Bw+lvyH  +  Cw+lvw+2  "  dw+l  ’ 


with  i  being  greater  than  or  equal  to  w+2  in  the  remaining  equations. 

1 

4.2  Detailed  Procedure 

Part  I  -  Before  Transformation 

Using  Equations  (3*/5)>  (3*6)  and  (3*1 6)  computation  continues  until 


for  some  j  =  o',  v0  j '+j_  *  1  and  vo,j/  < 
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To  establish  an  accurate  time  for  the  beginning  of  the  phase -change 
and  as  a  baseline  for  Part  II,  the  (j7+l)st  =  Jth  row  is  re-determined 
using  a  smaller  time  step  k7  given  by 


k7 


(  1  "  } 
0,  j  +1  0,0 


y 


or  equivalently 


r__ 

Vu 


1  -  U„  ./  ^ 

^ 
u. 


■o,j'+i  "  "o,  y 


The  time  of  the  start  of  transformation  is  simply 


t  =  j7k  +  k7  . 
m 


Part  II  -  Boundary  Begins  to  Move 

As  first  step  in  Part  II  the  boundary  must  be  located.  Since  the 
temperatures  at  time  J  +  1  are  yet  to  be  determined,  (3*15)  cannot  be 
used  directly.  We,  therefore,  approximate  Ay  by 


Ay 


j+i 


r 

2N 


’C-1*5 


+  2v. 


1,J 


-  *5v. 


2  ,  Jy 


+  (hpf 


J+l 


7 


(k.k) 


The  first  term  inside  the  brackets  may  be  recognized  as  the  three - 
point  Lagrange  interpolation  expanded  in  the  neighborhood  of  the  moving 
boundary  with  w  =  s  =  0.  Multiplying  the  brackets  by  l/2  ensures  the 
Ay  is  not  over  estimated.  In  the  absence  of  any  boundary  displacement 
at  time  J  the  initial  estimate  of  the  position  of  the  boundary  is 


0 

y 


J+l 


=  Ay 


J+l 


(4.5) 


where  the  superscript  designates  the  number  of  the  iteration  n,  taken 
to  be  zero  initially. 
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Further 


0 


J+l 


is  decomposed  into  its  whole  and  fractional  components. 


0 

y  J+l 


h  CwVi +  s Vj)  . 


(4.6) 


At  this  point  (5.10),  (5*5)  and  (5*1 6)  are  solved  for  all  v±  J+1< 

It  should  be  emphasized  that  temperatures  v.  T  are  not  replaced  by  those 

i,«j 

at  v  ,  ,  until  the  iteration  for  the  (J+l)  step  has  been  completed, 
l,  J+l 

Using  s°T  and  the  appropriate  v.  T  '  s  a  new  estimate  of  the 

boundary  position  is  obtained  from  (5. 15).  This  then  becomes  AyJ'j+1 

from  which  y1J+1  is  derived.  As  an  improved  estimate  of  the  boundary 

position  the  n  and  (nfl)  iterative  values  of  yj+]_  are  combined  and  a  new 

s _  ,  and  wT  ..  are  found, 

J+l  J  +  l 


0 

y  j+i 


j+i 


2 


+  s 


J+ 


(4.7) 


Then  the  system  (4.2)  is  again  solved  and  the  process  is  repeated  a 
fixed  number  of  times  or  until  some  test  is  satisfied.  We  have  chosen 
the  test  that  the  absolute  value  of  the  ratio  of  the  difference  between 
two  successive  values  of  y  to  their  average  be  less  than  some  constant  0 
and  write 


n  n-1 

y  -  y _ 

(yn  +  yn~x ')  /2 


0  > 


(4.8) 


where  0  <  0  «  1. 


For  an  approximation  to  y  in  the  succeeding  time  step  J  +  2,  a 
linear  change  of  the  boundary  with  time  is  assumed, 


0 

y 


J+2  =  yJ+l  +  ^J+l 


(4.9) 
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where  yT  and  Ay  represent  values  of  the  quantities  obtained  in  the 
final  iteration  of  the  (j-f-l)st  step.  The  above  procedure  commencing 
with  (4.6)  is  then  repeated.  Computation  continues  until  t  s  t  , 

intux 

where  t  is  the  duration  of  the  input  pulse, 
max 

5.  ERROR  ANALYSIS 

With  assurance  of  stability  for  independent  h's  and  k*s,  error 
analysis  is  one  of  the  prime  considerations  of  a  numerical  method. 

This  subject  is  currently  under  extensive  study,  and  most  discussions 

9 

found  in  the  literature  are  of  a  qualitative  nature.  Two  main  sources 
of  error  generally  encountered  are  present  here  and  are  (a)  the  trunca¬ 
tion  error  arising  from  the  approximation  of  the  analytical  equations 
by  finite  differences,  and  (b)  the  round-off  error  due  to  use  of  only  a 
finite  number  of  decimal  places  in  the  arithmetic  operations  and  due  to 
the  fact  that  the  iterative  solution  is  only  continued  until  there  is 
no  change  out  to  a  certain  decimal  place.  These  errors  are  oppositely 
influenced  by  the  interval  size  h.  Decreasing  h  decreases  the  trunca¬ 
tion  error  but  in  general  increases  the  round-off  error. 

An  example  will  best  serve  to  check  the  accuracy  of  the  present 

solution.  For  a  semi -infinite  slab  and  a  constant  heat  input  the  time 

required  for  the  surface  to  reach  the  phase -change  temperature  can  be 

12 

compared  to  the  exact  analytical  solution.  These  times  are  shown  in 
Table  I  using  physical  data  for  aluminum  undergoing  a  phase-change  at 
the  melting  temperature.  It  is  readily  seen  that  with  decreasing  h 
at  constant  r  the  approximate  time  that  the  surface  reaches  the  melting 
temperature  more  nearly  approaches  that  of  the  exact  solution.  The 
error  in  the  approximate  solution  is  reduced  by  the  square  of  the  factor 
by  which  h  is  decreased.  Thus  if  h  is  diminished  by  a  factor  of  8,  e.g., 
in  going  from  h  =  .01  to  h  =  .00125,  the  error  will  decline  by  a  factor 
of  64. 
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TABLE  I 


Time  of  the  start  of  melting  for  a  semi -inf inite  solid  as  determined 

*  /  \ 

by  (a)  the  exact  analytical  solution  and  (b)  the  approximate  solution 
of  this  paper  for  varying  h  and  constant  r  =  0.4.  Based  on  physical 
data  of  Al. 

Exact  Solution  Approximate  Solution 

time,  sec  time,  sec  h  $  Error 


.2082  x  10"5 

.2475  X  10"5 

.01 

19 

.2180  x  10"5 

.005 

4.7 

.2107  x  10'5 

.0025 

1.2 

.2088  x  10-5 

.00125 

•  3 

See  Reference  12. 

K  =  .480,  c  =  .214,  p  =2.70,  L  =  94.5,  Tq 
H  =  2.07  x  10 4,  B  =  1.0,  f(T )  =  1.0. 


20,  Tl  =  660, 
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In  Figure  1  the  results  of  Table  I  are  extended  into  Part  II  involv¬ 
ing  the  moving  boundary.  Here,  however,  it  is  no  longer  possible  to  cal¬ 
culate  an  exact  solution  and  the  approximate  results  must  be  compared  to 
the  steady  state  solution  (l)  given  by 


y  = 


E  (t  -  tj) 

+  c  (' T.  -  T  "V! 
V  L  o  j 


(5.1) 


where  in  addition  to  the  quantities  already  defined  t  is  the  duration 
of  the  pulse  and  t  is  the  time  required  to  reach  the  transformation 
temperature . 

Equation  (5.1)  gives  an  upper  limit  to  the  penetration  of  the 

-2 

boundary,  which  in  the  time  t  =  .24  x  10  sec  is  0.0725  cm.  Even  for 
the  least  favorable  choice  in  h  the  solution  is  reasonable.  Except  in 
the  very  short  time  after  transformation  has  begun,  the  boundary  motion 
is  essentially  linear  with  time. 

It  is  not  clear  from  the  preceding  discussion  how  one,  short  of 
trial  and  error,  goes  about  making  an  appropriate  selection  of  h.  Time 
and  accuracy  must  be  simultaneously  considered  and  some  optimum  found 
between  the  two.  Rather  than  select  h  first  it  is  more  meaningful  to 
choose  k  first  and  then  find  the  desired  h  by  h  =  (k/r)  7  .  This  is 
necessary  because  any  choice  in  h  no  matter  how  small  might  not  be 
applicable  due  to  the  short  duration  of  an  input  pulse.  For  a  nano¬ 
second  long  input  pulse  as  encountered  in  the  q-switched  laser  mode  an 
h  of  say  .00125  would  yield  a  k  of  0.625  x  10  for  r  =  0.4  which 
exceeds  the  entire  length  of  the  pulse.  In  choosing  a  k  one  may  be 
guided  by  a  rule  of  thumb  in  which  the  total  time  of  the  pulse  is 
divided  into  about  one  thousand  kfs  i.e.,  max  j  =  1000.  Then  a  back 
calculation  can  be  made  to  find  an  approximate  h  which  can  also  be 
divided  into  B  to  give  some  whole  number  M. 


In  implicit  methods  stability  is  ensured  for  all  r>  0*  [Ref,  93  0,3391, 
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When  the  input  energies  are  very  large  as  in  the  case  of  the  laser 
the  time  taken  for  the  material  to  reach  the  transformation  temperature 
is  much  shorter  (sometimes  by  a  factor  of  a  thousand  of  more)  than  the 
duration  of  the  entire  pulse.  It  is,  therefore,  necessary  to  use 
different  grid  sizes  for  Parts  I  and  II.  If  a  single  grid  size  is 
chosen  say  on  the  basis  of  the  high  resolution  required  in  Part  I  the 
computational  times  required  for  Part  II  would  become  prohibitive.  If 
on  the  other  hand,  the  choice  of  grid  sizes  is  made  on  the  basis  of  the 
total  pulse,  the  approximation  of  the  time  of  the  phase -change  will 
contain  a  ridiculously  large  error.  A  rough  idea  for  the  onset  time  of 
the  transformation  may  be  obtained  from  an  analytical  solution  of  the 
heat  equation.  On  this  basis  a  reasonable  selection  of  k  and  hence  h 
may  be  made  for  Part  I.  The  h  selected  in  this  manner  may  be  so  fine, 
however,  that  for  the  particular  thickness  of  solid  B,  M  will  become  an 
exceedingly  large  number.  Insofar  as  the  heat  flow  is  concerned,  only 
a  small  fraction  of  the  total  M  points  in  the  solid  may  register  any 
change  in  temperature  in  the  short  times  considered,  the  solid  behaving 
as  semi -inf inite  beyond  a  certain  depth.  'Thus  by  either  reducing  B  or 
equivalently  reducing  M  by  some  arbitrary  factor  the  times  for  calcula¬ 
tion  in  Part  I  do  not  become  prohibitive.  When  Part  II  is  entered  a 
much  larger  h  may  be  used.  To  avoid  the  possibility  of  error  the  final 
temperature  v^  may  be  printed  as  part  of  the  output  at  the  ends  of 
Parts  I  and  II.  If  v„  is  the  same  as  the  initial  temperature  of  the 
solid  then  the  slab  has  remained  semi -inf inite  throughout  the  length  of 
the  calculations. 


6.  AN  EXAMPLE  FOR  LASER -INDUCED  DAMAGE 

In  suggesting  that  the  present  problem  may  apply  to  laser -induced 
damage  in  materials  we  are  cognizant  of  the  pitfalls  involved.  The 
interaction  of  the  laser  with  the  target  is  a  complex  phenomenon. 
Several  elements  contribute  to  this  complexity.  It  is  not  certain 
whether  the  plasma  created  in  front  of  the  target  remains  transparent 
to  the  beam  throughout  the  duration  of  the  pulse.  It  may  be  a  safe 
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assumption  to  make  at  low  incident  power  densities  but  not  at  high 
incident  power  densities  such  as  those  delivered  by  the  laser  in  the 
q-switched  mode.  In  addition  to  the  plasma  the  target  emits  globules 
of  transformed  material  which  may  further  act  to  absorb  the  incoming 
beam.  Another  important  consideration  is  the  reflectance  of  the  beam 
by  the  target  which  may  be  a  strong  function  of  the  depth  of  hole,  i.e., 
in  the  initial  stages  of  heating  and  penetration  highly  reflecting 
materials  such  as  aluminum  may  absorb  only  a  small  portion  of  the 
incoming  radiation,  while  at  a  later  stage  or  as  the  hole  deepens,  the 
entire  incoming  beam  may  be  absorbed  by  being  trapped  in  the  hole.  The 
model,  furthermore,  makes  no  allowance  for  the  possibility  of  superheating 
or  for  the  co-presence  of  vapor  and  liquid.  Any  attempt,  however,  to 
construct  a  mathematical  model  of  laser  damage  which  takes  all  the  pre¬ 
ceding  factors  into  account  would  prove  an  enormously  difficult  task 
because  a  clear  understanding  of  their  relative  significance  is  currently 
lacking. 

In  Table  II  and  Figure  2,  data  are  shown  comparing  the  depth  of  hole 
predicted  to  that  actually  observed  in  A1  of  semi -infinite  thickness. 

Case  1  in  Table  II  is  for  a  constant  input  energy,  whereas  Case  2  if  for 
the  variable  pulse  shown  at  the  top  of  Figure  2.  The  lower  portion  of 
Figure  2  is  a  plot  of  the  position  of  the  boundary  with  time  for  this 
variable  pulse  with  the  maximum  y  being  the  entry  in  Table  II. 

In  cases  1  and  2  the  transformation  is  assumed  to  occur  at  the 
temperature  of  vaporization,  the  possibility  of  melting  being  totally 
ignored.  Ignoring  melting  when  vaporization  is  involved  is  not  a 
serious  error  since  the  heat  of  vaporization  and  melting  are  so  vastly 
different,  the  former  being  more  than  twenty-seven  fold  greater  than 
the  latter  in  Al. 

The  discrepancy,  assuming  full  absorption  of  the  energy  of  the 
beam  by  the  target,  is  quite  large.  A  second  calculation  has  been 
done  for  each  case  at  about  half  of  this  energy  on  the  basis  that  much 
of  the  incident  pulse  is  reflected.  Specimens  of  Al  treated  with  #300 
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TABLE  II 

DEPTH  OF  HOLE  PRODUCED  IN  A1  BY  LASER  PULSE 


Case 

Description  of 

Pulse 

Power  Density 

Duration 

Calculated 

Observed* 

joules/ cm^-sec 

sec 

cm 

cm 

1 

8.55  x  lo6 

.6  x  10  ^ 

1.45  x  10'1 

•78  x  IO"1 

4.16  x  106 

.6  x  10“5 

.695  x  10_1 

2 

1  x  109  (peak)4 

44  x  io "9 

6.05  x  10'4 

,  -4 

5.6  x  10 

•5  x  109 

44  x  IO-9 

2.87  x  10"^ 

* 


+ 


See  Reference  4. 

See  Figure  2  for  pulse  shape. 


2b 


y  X  I03 
(cm  ) 


FIG. 2.-B0UNDARY  MOTION  y  AS  A  FUNCTION  OF  THE  TIME  t  FOR 
INPUT  PULSE  Hf (t )  (SHOWN  ON  TOP  GRAPH).  MATERIAL 
IS  A I  AT  VAPOR  TEMPERATURE. 
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grit  emery  paper  have  been  found  in  this  laboratory  to  reflect  up  to  60 
per  cent  of  the  incident  light  at  room  temperature.  While  large 
reflectance  is  a  reasonable  assumption  during  the  heating-up,  it  is  not 
known  whether  it  will  remain  at  or  near  the  initial  level  for  the  entire 
length  of  the  pulse.  Better  agreement  with  the  experiment  results  when 
making  the  assumption  of  a  constant  reflectance  throughout  the  pulse 
duration. 


7.  CONCLUSIONS 

An  implicit  numerical  scheme  has  been  used  successfully  to  solve 
the  heat  conduction  problem  for  an  ablating  solid.  The  method  intro¬ 
duces  simplicity  and  accuracy  in  the  calculations  and  allows  for 
economic  use  of  computer  time.  Further,  the  solution  is  valid  for  a 
variable  input  energy  and  sample  thickness,  features  which  make  it 
especially  suitable  in  applications  related  to  the  laser  damage  of 
materials . 

In  the  example  given  for  laser-induced  damage  to  Aluminum  we  have 
purposely  avoided  manipulating  our  input  energies  etc.,  so  as  to  obtain 
a  better  agreement  with  experimental  observations.  Published  data  of 
laser  damage  to  targets  is  scarce  and  incomplete.  Information  concern¬ 
ing  surface  preparation,  the  spread  in  the  measured  input  energies  and 
depths,  and  the  shape  of  the  holes  produced  has  never  been  set  down  for 
the  same  samples.  In  certain  cases  one  cannot  be  certain  whether  a 
semi -inf inite  slab  was  used  as  reported  or  if  the  conditions  for  one- 
dimensional  heat  flow  were  truly  satisfied.  From  the  limited  experi¬ 
mental  observations  considered,  it  is  difficult  to  determine  how  widely 
the  proposed  model  can  be  applied  to  predicting  laser  damage  to  targets. 
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